# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import warnings
import numpy as np
import sympy as sm
from hysop.tools.htypes import first_not_None, check_instance
from hysop.symbolic import space_symbols
from hysop.symbolic.field import div, grad
from hysop.tools.warning import HysopWarning
from hysop.tools.sympy_utils import get_derivative_variables
[docs]
def collect_direction(expr, var):
is_tensor = isinstance(expr, np.ndarray)
if is_tensor:
R = np.empty_like(expr)
for idx in np.ndindex(*R.shape):
R[idx] = collect_direction(expr[idx], var)
return R
else:
args = ()
if isinstance(expr, sm.Symbol):
return expr
elif isinstance(expr, sm.Derivative):
if set(get_derivative_variables(expr)) != {var}:
return 0
else:
return expr
elif isinstance(expr, (sm.Number, int, float)):
return expr
elif isinstance(expr, sm.Expr):
for e in expr.args:
de = collect_direction(e, var)
args += (de,)
return expr.func(*args)
else:
msg = (
f"Unknown expression type {type(expr).__mro__}.\nExpression was: {expr}"
)
raise TypeError(msg)
[docs]
def split(F, fixed_residue, force_residue, coords):
from hysop.symbolic.relational import Assignment
check_instance(F, sm.Basic)
if isinstance(F, Assignment):
msg = f"Expression cannot be of type Assignment: {F}"
raise TypeError(msg)
res = {}
coords = first_not_None(coords, space_symbols)
for i, xi in enumerate(coords):
res[i] = collect_direction(F, xi)
try:
if force_residue is not None:
residue = force_residue
else:
residue = (F - sum(res.values())).simplify()
except AttributeError:
msg = "Failed to compute residue for expression {}."
msg = msg.format(F)
raise RuntimeError(msg)
count = sum((r != 0) for r in res.values())
if count > 0:
for i in res.keys():
res[i] += residue / count
if fixed_residue is not None:
for i, r in res.items():
if r != 0:
res[i] += fixed_residue
return res
[docs]
def split_assignement(expr, fixed_residue, force_residue, coords):
lhs, rhs = expr.args
res = split(rhs, fixed_residue, force_residue, coords=coords)
for i, ei in res.items():
if ei != 0:
res[i] = expr.func(lhs, ei)
return res
if __name__ == "__main__":
from hysop import Field, Box
from hysop.tools.sympy_utils import enable_pretty_printing, greak
from hysop.parameters.scalar_parameter import ScalarParameter
enable_pretty_printing()
box = Box(length=(1,) * 3)
frame = box.frame
U = Field(domain=box, name="U", is_vector=True)
W = Field(domain=box, name="W", is_vector=True)
alpha = ScalarParameter(name=greak[0]).s
u = U.s(*frame.coords)
w = W.s(*frame.coords)
expr = (alpha * grad(u, frame) + (1 - alpha) * grad(u, frame).T).dot(w)
# expr = div(np.outer(u,w), frame)
print(expr)
for i, x in split(expr, frame.coords).items():
print(i, x.simplify())